#Dissection simulations 
library(HeadR)  # loads data.table too.
library(SQUAREM)
library(lfe)
library(doParallel)
library(doRNG)
library(xtable)
#if you opened the session by doubleclicking the R project in the BLP folder then switch directories.
setwd("~/Dropbox/BLP_REStat/Replication/Rfiles")
setDTthreads(1)
rm(list=ls())
source("CF_MMX_functions.R")
source("CF_MMX_oly_functions.R")
source("CF_MMX_MLOG_functions.R")
library(Rcpp)
Rcpp::sourceCpp("../Cppfiles/crossDelta.cpp")
#
cl <- detectCores()-1
registerDoParallel(cores=cl)
getDoParWorkers()
MAXIT <- 1000
#
dol.units <- 1000 # count in 000s of dollars (for income variables)
#
#------------------------------------------------------------#
#Parameters common to all====
#------------------------------------------------------------#
seedn <- 777 #666
n.hh <- 1000 # households
ntb <- 0 # this is the ntb part of tau, not needed anymore 
tar <- 0.1 # this is the tariff: 10%
#dist.ave.pars <- c(0.2,0.2) # if random
dist.ave.pars <- c(0.1,0.1)
n.countries <- 3
n.models <- 90
models.perfirm <- 10
n.x <- 4  # dimensions of quality (number of measured model attributes) 
U0 <- 1  # outside good, in this case there will be constant x0 =1, 
#U0 <- 0 # no outside good, x1..x{n.x}
GINI <- FALSE
sigma.ly.factor <- 1
#n.reps <- 20
n.reps.target <- 1000
n.reps <- n.reps.target + 60
nv <- 3 #nrow(settings.var)
#
#------------------------------------------#
#  OG90 simulations ====
#------------------------------------------#
#reset the settings
OGpct <- 90
source("CF_MMX_settings_postcalib.R") # loads the relevant settings (including OGpct option)
#
#MCES
#
#Panel B (will be the new "(c)" in tables): MC approx of  Multiproduct Oly 
mix.var <- "MCES/MC"
n.firms <- n.models/models.perfirm
settings.var <- b.settings
pan.var <- "b"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
monte.func.mm(1,settings=settings.var,v=1,market=1) # test one rep
#create a dataframe with all four versions
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mm(parallel=T,settings=settings.var,v,market=1)}))
#format output
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#-------------------------------------------#
# Signal if there were convergence problems #
#-------------------------------------------#
files <- list.files(path="../Tables/CF_MM/MCES/MC", pattern="monte_allv", all.files=T, full.names=T)
for (f in files) {
  temp <- readRDS(file=f)
  cat(f,": ", sum(temp$converge.prob)," convergence failures, out of ",length(temp$converge.prob), "reps \n")
}
#
# MLOG
#
#Panel A: MC
mix.var <- "MLOG/MC"
n.firms <- n.models
settings.var <- a.settings.mlog
pan.var <- "a"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
monte.func.mlog.mm(1,settings=settings.var,v=2,market=1) # test one rep
#create a dataframe with all four versions
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mlog.mm(parallel=T,settings=settings.var,v,market=1)}))
#format output
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#Panel B: Multiproduct Oly 
n.firms <- n.models/models.perfirm
settings.var <- b.settings.mlog
pan.var <- "b"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
monte.func.mlog.mm(1,settings=settings.var,v=3,market=1) # test one rep
#create a dataframe with all  nv settings
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mlog.mm(parallel=T,settings=settings.var,v,market=1)}))
#format output
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#-------------------------------------------#
# Signal if there were convergence problems #
#-------------------------------------------#
files <- list.files(path="../Tables/CF_MM/MLOG/MC", pattern="monte_allv", all.files=T, full.names=T)
for (f in files) {
  temp <- readRDS(file=f)
  cat(f,": ", sum(temp$converge.prob)," convergence failures, out of ",length(temp$converge.prob), "reps \n")
}


#------------------------------------------#
#  OG10 simulations ====
#------------------------------------------#
#reset the settings
OGpct <- 10
source("CF_MMX_settings_postcalib.R") # loads the relevant settings (including OGpct option)
#
#MCES
#
#Panel B (will be the new "(c)" in tables): MC approx of  Multiproduct Oly 
mix.var <- "MCES/MC"
n.firms <- n.models/models.perfirm
settings.var <- b.settings
pan.var <- "b"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
monte.func.mm(1,settings=settings.var,v=3,market=1) # test one rep
#create a dataframe with all four versions
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mm(parallel=T,settings=settings.var,v,market=1)}))
monte.allv[,.(mean(change),mean(change.ces)),by=v] #quick check of results
#format output
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#Panel C: OLY approx of  Multiproduct Oly 
mix.var <- "MCES/OLY"
n.firms <- n.models/models.perfirm
settings.var <- c.settings
pan.var <- "c"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
#monte.func.cf.oly.mm(1,settings=settings.var,v=1,market=1) # this should always replicate the true change
monte.func.cf.oly.mm(1,settings=settings.var,v=3,market=1)
#create a dataframe with all four versions
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.cf.oly.mm(parallel=T,settings=settings.var,v,market=1)}))
#output results
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#-------------------------------------------#
# Signal if there were convergence problems #
#-------------------------------------------#
files <- list.files(path="../Tables/CF_MM/MCES/MC", pattern="monte_allv", all.files=T, full.names=T)
for (f in files) {
  temp <- readRDS(file=f)
  cat(f,": ", sum(temp$converge.prob)," convergence failures, out of ",length(temp$converge.prob), "reps \n")
}

files <- list.files(path="../Tables/CF_MM/MCES/OLY", pattern="monte_allv", all.files=T, full.names=T)
for (f in files) {
  temp <- readRDS(file=f)
  cat(f,": ", sum(temp$converge.prob)," convergence failures, out of ",length(temp$converge.prob), "reps \n")
}
#
#MLOG
#
#Panel B: Multiproduct Oly 
mix.var <- "MLOG/MC"
n.firms <- n.models/models.perfirm
settings.var <- b.settings.mlog
pan.var <- "b"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
monte.func.mlog.mm(1,settings=settings.var,v=3,market=1) # test one rep
#create a dataframe with all  nv settings
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mlog.mm(parallel=T,settings=settings.var,v,market=1)}))
#format output
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#-------------------------------------------#
# Signal if there were convergence problems #
#-------------------------------------------#
files <- list.files(path="../Tables/CF_MM/MLOG/MC", pattern="monte_allv", all.files=T, full.names=T)
for (f in files) {
  temp <- readRDS(file=f)
  cat(f,": ", sum(temp$converge.prob)," convergence failures, out of ",length(temp$converge.prob), "reps \n")
}


#------------------------------------------#
#  OG50 simulations ====
#------------------------------------------#
#reset the settings
OGpct <- 50
source("CF_MMX_settings_postcalib.R") # loads the relevant settings (including OGpct option)
#
#MLOG
#
#Panel B: Multiproduct Oly 
mix.var <- "MLOG/MC"
n.firms <- n.models/models.perfirm
settings.var <- b.settings.mlog
pan.var <- "b"
#
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
tau.grid.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do%  make.tradecosts(n.countries,n.firms,n.models,v,settings=settings.var))
tau.grid.allv[,tariff.change:= tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
set.seed(seedn)
monte.func.mlog.mm(1,settings=settings.var,v=3,market=1) # test one rep
#create a dataframe with all  nv settings
#set.seed(seedn)
system.time(monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mlog.mm(parallel=T,settings=settings.var,v,market=1)}))
#format output
monte.output.mm(simdata = monte.allv,mix = mix.var,pan=pan.var)
rm(monte.allv,tau.grid.allv,co.own.mat) # for security
#
#-------------------------------------------#
# Signal if there were convergence problems #
#-------------------------------------------#
files <- list.files(path="../Tables/CF_MM/MLOG/MC", pattern="monte_allv", all.files=T, full.names=T)
for (f in files) {
  temp <- readRDS(file=f)
  cat(f,": ", sum(temp$converge.prob)," convergence failures, out of ",length(temp$converge.prob), "reps \n")
}


